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Abstract. We investigate the Lagrangian perturbation theory of a homogeneous and 
isotropic universe in the non-relativistic limit, and derive the solutions up to the fourth 
order. These solutions are needed for example for the next-to- leading order correction of the 
(resummed) Lagrangian matter bispectrum, which we study in an accompanying paper. We 
_i- | focus on flat cosmologies with a vanishing cosmological constant, and provide an in-depth 

^J ■ description of two complementary approaches used in the current literature. Both approaches 

are solved with two different sets of initial conditions — both appropriate for modelling the 
large-scale structure. Afterwards we consider only the fastest growing mode solution, which 
is not affected by either of these choices of initial conditions. Under the reasonable approxi- 
mation that the linear density contrast is evaluated at the initial Lagrangian position of the 
S^ | fluid particle, we obtain the nth-order displacement field in the so-called initial position limit: 

the nth order displacement field consists of 3(n-l) integrals over n linear density contrasts, 
and obeys self-similarity. Then, we find exact relations between the series in Lagrangian and 
Eulerian perturbation theory, leading to identical predictions for the density contrast and 
the peculiar- velocity divergence up to the fourth order. 
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1 Introduction 

In the last years several analytic techniques have been proposed in order to study the inho- 
mogeneities of the large scale structure (LSS) [1-7]. The basic idea of them is to solve the 
equations for an irrotational and pressureless fluid of cold dark matter particles in terms of 
a perturbative expansion. In the standard scenario the density and velocity field of the fluid 
particle are the perturbed quantities. Thus the validity of the perturbative series depends 
on the smallness of these fields. This approach is called Eulerian (or standard) perturbation 
theory (SPT), since the equations are evaluated as a function of Eulerian coordinates [8]. 
Subsequent gravitational collapse leads to highly non-linear structures in the universe like 
galaxies, clusters of galaxies, etc., i.e., regions where the local density field departs signifi- 
cantly from the mean density. As a result, the series in SPT breaks down. This situation was 
already realised in 1969 by Zel'dovich, who proposed an approximate solution which is above 
all applicable to the highly non-linear regime by a Lagrangian extrapolation of the Eulerian 
linear solution, inspired by the exact solution for inertial systems [9-12]: the Zel'dovich ap- 
proximation (ZA). In general, the ZA can be derived from the full system of gravitational 
equations and forms a subclass of solutions of the Lagrangian theory of gravitational insta- 
bility (i.e., Lagrangian perturbation theory; LPT) [3, 13-16]. In LPT the only perturbed 
quantity is the gravitational induced deviation of the particle trajectory field from the homo- 
geneous background expansion. Stated in another way, LPT does not rely on the smallness 
of the density and velocity fields, but on the smallness of the deviation of the trajectory 
field, in a coordinate system that moves with the fluid. It can be shown that this implies 
a weaker constraint on the validity of the series and hence can be maintained substantially 
longer during the gravitational evolution (see the thorough discussion in [17-19]). Addi- 
tionally, to obtain an SPT series one basically has to approximate the continuity- and the 
Euler-equation order by order, so that, strictly speaking, mass- and momentum-conservation 
are not fulfilled. In LPT, on the other hand, the Jacobian of the transformation from the 
Eulerian to the Lagrangian frame is approximated and so is the precise localisation of the 
fluid element, whereas the continuity- and Euler-equation are still exactly solved. 

General perturbation and solution schemes to any order of Lagrangian perturbations 
on any FLRW background have been given in the review [20], based on explicit evalua- 
tions of the general first-order scheme including rotational flows [3, 13, 21], and the general 
second-order scheme for irrotational flows [15] . The third-order scheme for irrotational flows 
with slaved initial conditions (i.e., for an assumed initial parallelism of peculiar- velocity and 
peculiar-acceleration) is given in [16], and the fourth-order scheme for this subclass of the 
general solution has been derived in [22] (see further below for an explanation). Lagrangian 
perturbation theory has also been extended to include pressure [23], and the series can be 
derived from exact integrals of longitudinal and transverse parts [24, 25]. Extensive compar- 
isons of LPT results against N-body simulations can be found in [2, 19, 26-32]. 



In the present paper, we explicitly reexamine the Lagrangian framework in two different 
representations and evaluate them up to the fourth order. For both representations we choose 
a different set of initial conditions, which can be labeled as 'Zel'dovich type', since only one 
initial potential has to be chosen instead of two in the general case. Note that the fastest 
growing mode solution is not affected by either of these choices. At this point the reader 
may ask, what is the point in deriving higher-order solutions for the purpose of modelling the 
LSS. First of all, the fourth-order solution is needed for the next-to- leading order correction 
to the LPT matter bispectrum, which we calculate in an accompanying paper. From a 
theoretical point of view one also expects a match between SPT and LPT under certain 
circumstances. In reference [5] the equivalence of SPT and LPT is shown, if one sums up 
the perturbative solutions up to the third order. However it is not clear a posteriori if this 
matching between SPT and LPT occurs at the fourth-order as well, because the convergence 
of the LPT and SPT series need not to behave equally. Furthermore, there is a growing 
interest to apply higher-order solutions in the context of resummation approaches, and we 
consider an explicit demonstration and a thorough comparison with the SPT series useful. 
Also, note that resummation techniques in the LPT framework are directly feasible rather 
than complex scenarios in SPT [33]. 

Although the subject of this paper is quite technical, we try to keep it as readable 
as possible, e.g., we restrict our calculations to an Einstein-de Sitter (EdS) universe. The 
organisation of this paper is as follows: in section 2 we derive the evolution equations step by 
step and confront two complementary ways of how to deal explicitly with LPT calculations. 
We do so to shed light on two different looking formalisms used in the current literature. 
Then, in section 3 we mention and explain our choice of initial conditions. In section 4 and 5 
we show the results in both formalisms. Afterwards, we prepare our solutions to be used in 
Fourier space and derive relations between the SPT and LPT series in section 6. Finally, in 
section 7 we give a discussion and conclude. Our notation is introduced and defined in the 
text, but is also summarised in table 1. 

2 Systems of equations 

According to the ACDM model and its current success in treating most of the problems in 
observational cosmology [34], we live in a statistically homogeneous and isotropic universe. 
The universe is expanding, thus the mean density ~p(t) is diluting. But this global effect 
cannot compete with the gravitational potential locally. Hence local density fluctuations 
5(r,t) are the source of gravitational collapse, which leads to the observed LSS. We can 
define the above quantities in terms of the full density p(r,t) as 

p(r,t)=p(t)[l + 5(r,t)}, (2.1) 

where p is given by the assumed homogeneous background density in a Newtonian model 
with hypertorus topology [35]. In this chapter we set up the evolution equations that are 
sourced by 

• 2.n.l: the full density p(r,t), which we label with "full-system approach", 

• 2.n.2: and for the density contrast 5(r,t), the "peculiar-system approach", 

step by step (n = 1, • • • , 4) and independent from each other. Readers who are only 
interested in the final equations may go directly to page 5: Eqs. (2.21-2.23) are the equations 



in the full-system approach, whereas Eqs. (2.29-2.31) refer to the peculiar-system approach. 
The perturbation equations are shown in section 2.4. Finally, in section 2.5, we clarify errors 
and common misunderstandings in the current literature. 

2.1 Eulerian equations 

2.1.1 Full-system approach 

Let us briefly go through the derivation of the equations of motion in the Lagrangian de- 
scription. For simplicity we focus on flat cosmologies with a vanishing cosmological constant 
although a more general implementation is straightforward [3, 13, 36]. As usual we denote the 
density, the velocity and the acceleration fields by p, v, and g, respectively. In a non-rotating 
(Eulerian) frame with coordinate r at cosmic time t, the equations for self-gravitating and 
irrotational dust are [37]: 

^A + V r -[p(r,t)v(r,t)]=0, (2.2) 

e ijk d rj v k (r,t) = 0, (2.3) 

e ijk d rj g k {r,t) = 0, V ' (2.4) 

V r -v(r,t) = -4irGp(r,t), with dv/dt = v = g and p > , (2.5) 

where Einstein summation over the spatial (Eulerian) components is implied. Eq. (2.2) is the 
continuity equation and denotes mass conservation, Eq. (2.3) states the irrotationality of the 
velocity, and should be viewed as an additional constraint to the field equation that requires 
an irrotational acceleration field, i.e., to Eq. (2.4). The divergence of the field strength, here 
with Euler's equation inserted, is linked to the density source in Eq. (2.5). Note that we make 
use of the convective time derivative, i.e., d/dt := d/dt\ r + v ■ V r , which we denote by an 
over dot. 

2.1.2 Peculiar-system approach 

Eqs. (2.2)-(2.5) are written in terms of the full density p(r, t), thus including the homogeneous 
and isotropic deformation of an expanding universe p(t) (= p /a s ) for a matter dominated 
universe). However, it is also possible to construct a set of equations [1, 13, 38] where Poisson's 
equation is only sourced by the density contrast S, which is linked to the full density in the 
following way: p(x,t) = ~p(t) [1 + 5(x,t)]. x denotes the comoving distance and is related to 
the physical distance as r = ax, where a(t) is the cosmic scale factor. The Poisson equation 
then reads [13, 39-41]: 

A x $(x,r t ) = a(r t )8(x, V ), a = 6/(r, 2 + k), (2.6) 

where we have switched to super confer mal time 7] = \J — k(l— f])" 1 ' 2 (we denote its derivatives 
with d/d?] = ') [38]. For an EdS universe we have the simplification a?drj = d£, with 
a = (rjo/rj) 2 , and a = 6/rj 2 . The peculiar-evolution equations can then be written as (compare 
with Eqs. (2.2-2.5)): 

^^ + V x ■ {[1 + S(x, V )} u(x, V )a} = 0, (2.7) 

£ijkd Xj u k (x,r]) = 0, 
e ijk d Xj u' k (x,rj) = 0, 
v * • 9 P cc( x i 7 l) = -a(T]) 6(x,rj) , with du/drj = u' = g pec 



(i = 1,2,3) 


(2.8) 
(2.9) 


and 5 > — 1 , 


(2.10) 



Eq. (2.8) states the irrotationality of the peculiar- velocity u = adx/dt, Eq. (2.9) the irro- 
tationality of the (rescaled) peculiar-acceleration g = d 2 x/d 2 r], and Eq. (2.10) links the 
acceleration to the density contrast field. 

2.2 From the Eulerian to the Lagrangian framework 

The two set of equations, namely Eqs. (2.3-2.5) and Eqs. (2.8-2.10) are equivalent but have 
technical subtleties with their pros and cons, which we point out in the following Lagrangian 
description. We now briefly recall the corresponding Lagrangian systems that have been 
introduced in [37] and [13] in the full-system approach, and in [13], appendix A, for the 
peculiar-system. Note that the Lagrangian description can be formulated in the same La- 
grangian coordinates in both systems, which is the reason why the peculiar-system is es- 
sentially redundant. We nevertheless have chosen to confront the two approaches for the 
purpose of assisting work that deals with either of the two representations. Additionally, the 
peculiar-system approach is useful in order to link the LPT series to its counterpart in SPT, 
which we do so in section 6. 

2.2.1 Full-system approach 

It is useful to transform from Eulerian coordinates r,; to Lagrangian coordinates qi (i = 
1,2,3), and we start with the transformation of Eqs. (2.2-2.5). We introduce integral curves 
r = f(q, t) of the velocity field 

%-v. (2.11) 

and set the initial position at time to to 

f(q,t )=:q. (2.12) 

The Jacobian of the transformation can be written as 

Jij(q, t) := fi d , J := det[J^] , (2.13) 

where commas denote partial derivatives with respect to Lagrangian coordinates q. The 
formal requirement J > guarantees the existence of regular solutions and the mathematical 
equivalence to the Eulerian system [20]. The continuity equation, Eq. (2.2), is then integrated 
to yield p = p(q, to)/ J, and with the usage 1 of g = f we cast Eqs. (2.3-2.5) into [37] 

£i jk Jrj l Jkl = 0, (2.14) 

SijkJ^Jk^O, (2.15) 

Jf/Jji = -^GpoJ- 1 , p = p(q, t ) . (2.16) 

Note that Eq. (2.15) follows directly from Eq. (2.14). Hence from the irrotationality of the 
velocity field we can conclude the irrotationality of the acceleration field (but not vice 
versa!). This is shown in [15] and is valid for any integral curves, and for the perturba- 
tive solutions at each order as well. With the inverse Jacobian, i.e., J^- = l/Jadj[Jjj] = 



1 We assume the equivalence of the acceleration field and the gravitational field strength, i.e., Einstein's 
equivalence principle. 



1/ (2J) Eu m Ej pq JpiJ qm , and the use of Eq. (2.13) we have: 



Ji,n £njkjl,jjl,k = U j 
Ji,n £njkjl,jjl,k = U, 
£ilm£jpqjj,ijp,ljq,m — o7T(_rpo • 

To solve Eqs. (2.17-2.19) we impose the following ansatz for the trajectory 

f(q, t) = a(t) q + p(q, t) =>- f i;j = a5ij + p itj , 



(2.17) 
(2.18) 
(2.19) 



(2.20) 



where aq stands for the homogenous-isotropic background deformation, and p is the perturba- 
tion — induced by gravitational interaction. Plugging Eq. (2.20) into Eqs. (2.17-2.19) we finally 
obtain: 



(a 



d 

d7 
d_ 2 
dt 2 

2 d_ 2 
di 2 



aa) Eijk Pk,j = a e^k pij pi^ + Pi,n £njk 

~-njk 



Pl,j Pl,k ~ (a 



d 

d7 
d 2 



a ) Pk,j 



PijPi,k-(a^2-a)Pk,j 



(a 2 -772 ~ da ) £ ijk Pk,j = a e ijk p t)j p l)k + p i>n e r 

{a 2 ^72 + 2da ) Phi + a (Pi,i PjJ ~ Phi Pj,i) + 2 &'* Pi'i ~ Pi 'i p ^ + p tj &>* 

= -(47rGp + 3aa 2 



(2.21) 
(2.22) 

(2.23) 



In the above equations we have defined the co-factor element £>f~ = l/2£iim£jpqPp,iPq,m- 
Eqs. (2.21-2.23) with po > form a closed set of Lagrangian evolution equations in the full- 
system approach. 

2.2.2 Peculiar-system approach 

Analogous considerations lead to a similar set for Eqs. (2.8-2.10): the comoving trajectory 
field is x = F(q,rj), and with g = F" it follows that: 



»pcc 



^■ilm^jpq^ j i"p,l^ q,m = aQOJ , 



(2.24) 
(2.25) 
(2.26) 



where J = det [i^j], and, as before, a prime denotes a derivative with respect to supercon- 
formal time. For the set of Eqs. (2.24-2.26) we impose the ansatz: 



F(q,r]) =q + ^(q,rj) => F, 



> -j 



Sij + ty , 



(2.27) 



where the crucial difference in Eq. (2.27) with respect to Eq. (2.20) is the missing factor of 
a, and F = f /a links the comoving to the physical trajectory field. In the peculiar-system 
approach it is common to choose mass conservation by the following constraint: 



p(x, rj) d 3 x = p(rj) d 3 q , 



1/J* 



1. 



(2.28) 



Explicitly, in doing so we either restrict ourselves to a specific class of initial conditions [3] 
or assumed initial quasi-homogeneity. We shall discuss this issue later in section 5. 



Plugging Eq. (2.27) into Eqs. (2.24-2.26) we have: 

£ ijkK,j = e»ifc*U*U + **.« £ ^jk (*W*U " *fcj) . ( 2 - 2 9) 

£ iifc*fcj = £yfc*U*i!fc + *«.« e »i* (*Ij*£* " *kj) < (2- 3 °) 

0(77)^ - 1] = [(1 + \P M )<% - *ij + ^] *£ • (2.31) 

Similar to above we have defined \P? • = 1/2 eii m £jpq^ Pi i \P 9)m . Eqs. (2.29-2.31) with J F > 
is the closed set of Lagrangian evolution equations in the peculiar-system approach. 

2.3 The perturbation ansatz 

The full- and peculiar-systems are highly non-linear, thus it is common to seek approximate 
solutions in terms of perturbative series. 

2.3.1 Full-system approach 

First we proceed with the ansatz for the full-system, i.e., Eqs. (2.21-2.23). We expand the 
perturbation p(q, t) into a series, and factorise out the spatial and temporal dependence. 
Additionally, we decompose the inhomogeneous deformation p(q, t) of the nth order into 
purely longitudinal contributions, which we denote by ^f^ n '(q) = W q (j)^ n ' {q) , and purely 
transverse contributions, 2 denoted by T^ n '{q) = V q x A^ n '(q); the temporal parts of the 
nth-order correction are denoted by q n (t) or q^(t): 

p(q, t) = eqy{t) *«(q) + e 2 q 2 (t) *( 2 )(q) + e i q 3 (t) *( 3 )(g) + e%(t) *^ (q) 

+ e i q nt)T^(q) + e\J(t)T^(q). (2.32) 

The parameter e is supposed to be small and dimensionless. In the most general cases, 3 
nontrivial solutions of the irrotationality condition (i.e., V q x T^ n ' 7^ 0) of Eq. (2.17) and 
Eq. (2.18) occur the first time at the order of e 2 and are henceforth required for higher-order 
solutions. It should be pointed out that in addition to longitudinal (i.e., potential) modes, 
Lagrangian transverse modes also affect the growth of density perturbations. Additional to 
the above perturbation ansatz, we set the initial density po (on the RHS in Eq. (2.19)) to be 

P (q,to) = p(to)[l + e8(q,to)} , (2.33) 

without loss of generality, where p(io) an d S(q,to) denote the initial background density and 
the initial density contrast, respectively. In section 4 we shall need the above equation in 
order to set up the initial data. 

2.3.2 Peculiar-system approach 

The perturbative treatment of \1/ is quite analogous to the above. However, because of 
the LHS in Eq. (2.31), we need the explicit expansion of the Jacobian as well. This is 
clearly a (technical) disadvantage in comparison with the evolution equations in the full- 
system approach, where the Jacobian cancels out (but only in the absence of a cosmological 
constant). 



2 The Lagrangian transverse parts are mandatory in the Lagrangian frame to guarantee irrotationality in 
the Eulerian frame. 

3c 



See the general second-order solution for irrotational flows given in [15]. 



The Jacobian of the transformation between the Eulerian and Lagrangian frame depends 
on the displacement field \&(g, 77), 

J F (q, n) = det [6ij + Vij] = 1 + %,i + - [Vifljj - %,j^j,i] + det [tfy] . (2.34) 

Note that this is an exact relation, i.e., valid for the exact displacement field SI/. As stated 
above, the exact *&(q, rj) is expanded in a series. As we shall see, the spatial parts of the per- 
turbations agree with the spatial parts in Eq. (2.32) at each order, so we keep the previously 
introduced notation for them, but relabel the temporal parts to D(rj), E(rj), etc.: 

9(q,n) = eD{ V ) ^ { - 1 \q)+e 2 E{r ] )^ { - 2 \q) + e 3 F(T])^ 3 \q) + e 4 G(T])^ 4 \q) 

+ e 3 F T ( V )T^(q)+e 4 G T (7 1 )T^(q).(2.35) 

In the following we suppress the explicit temporal and spatial dependences. From Eq. (2.35) 
we can approximate the Jacobian by: 

jF^ v ) = l + eD^ +e 2 D 2 ^ ] +e 2 E^ ] + e 3 F^ + 2e 3 DE^ 2) 
+e 3 D 3 ^ + 2e 4 DF^ 3) - e 4 DF T ^T^ 

+e 4 G^ + e 4 E 2 /4 2 > + \e 4 D 2 Ee ikl e jmn ^ m *«*$ , (2 . 36) 

and the ath scalars fia = (ia (q) are defined by 

tf> = *if , (2.37) 

fi^U^M}-*®*®), (2-38) 



2 V l ' 1 J > J l >3 ■?>* 



,W 



>) 



1 



(n) m (n) , T ,(n) 



^ ee det [*£| = - e M e jmn V^ZKi , (2.39) 

and specifically 

^Kn) _ 1 ^M^(n) _ ^(n)) . (2.40) 

Note that [r™ = /^ for any tensor tyf 1 - and not only for longitudinal fields as pointed 
out incorrectly in [38], since it only consists of interchangable dummy indices. 

The above scalars will also be used for the spatial parts of the nth-order perturbations 
in p, since they are identical. 

2.4 The perturbation equations in Lagrangian form for an EdS background 

As before, we start with the perturbation equations in the full-system approach. We first 
concentrate on solving the source equation, Eq. (2.23), and write down the perturbative 
irrotationality condition for the velocity field, Eq. (2.21). Afterwards, in section 2.4.2, we 
perform the same steps for the peculiar-system. 



2.4.1 Full-system approach 

Inserting the ansatz, Eq. (2.32) and Eq. (2.33) into Eq. (2.23), we obtain the following set of 
equations to be solved: 

e°{4vrGp + 3aa 2 } =0, (2.41) 

e 1 | [a 2 91 + 2da qi ] ^ + 4ttGp 5 } = , (2.42) 

e 2 | [a 2 q 2 + 2daq 2 ] nf ] + [2qiqia + dqf] $ } = , (2.43) 

e 3 | [a 2 q 3 + 2daq 3 ] pf ] + 2 [d qi q 2 + a qi q 2 + aq\q 2 ] i$' 2) + 2>q\ql^ ] } = , (2.44) 

2 

e 4 {[a 2 g 4 + 2aa 94 ] ^ + [dqj + 2o&ffi] A^ + e ilm e jpq [g 2 y + §19192] *$ *$**# 
+ 2 [agi93 + «9l93 + a?i9 3 ] $® ~ [aqiQa + aq\qj + 091 9?] *$ T^ } } = . (2.45) 

Note that in the fourth-order part there is the occurrence of the transverse perturbation T^ ' ; 
thus, in order to solve this equation we have to constrain it with the irrotationality condition. 
Inserting the ansatz into Eq. (2.21), the third- and fourth-order parts are: 

e 3 { [aql - aqj] e ijk T$ - [ qi q 2 - q lQ2 ] £i ifc ¥$*g } = , (2.46) 

e 4 { [aql - aqj] e ijk T$ - [9193 - 9i9s] %^f,^g 

- [9i9 3 T - QiqI] eijMpfi! } = ° ■ ( 2 - 47 ) 

As mentioned above, the irrotationality of the acceleration field follows from the irrotation- 
ality of the velocity field, as can be easily checked by time differentiating Eq. (2.46) and 
Eq. (2.47) and comparing it directly with the perturbation equation of Eq. (2.22). 

2.4.2 Peculiar-system approach 

The peculiar equations, by construction, start at the order 0{e 1 ). Inserting the ansatz, 
Eq. (2.35), and the Jacobian, Eq. (2.36), into the source equation Eq. (2.31), delivers: 

e 1 {[D"-aD]/Jp}=0, (2.48) 

e 2 { [E" - aE] ^ - [aD 2 - 2DD"] $ } = , (2.49) 

e 3 { [F" - aF] (j® + 2 [D"E + DE" - aDE] (j$> 2) 

+ [3D"D 2 - aD 3 ] 4 1} } = , (2.50) 
e 4 \[G" - aG] ^ + [2E"E - aE 2 ] ^ 



"T £ilm,£jpq 



D 2 E" „ D 2 E 
+ D"DE - a 



^W^(i) $( 2 ) 

p,l q,m *j,i 



+ 2 [D"F + DF" + aDF] $' 3) - [D"F T + DF% - aDF T ] 9^ljf\ = . (2.51) 



With the irrotationality condition of the peculiar- velocity, Eq. (2.29), we obtain: 

e 3 {F^s ijk TQ - [DE> - D'E] ey fc *«*{5} = 0, (2.52) 

^ {g'tS^T^ - [DF> - D'F] e ***g*g 

- [DFlr - D'F T ] ey^T® } = . (2.53) 

Again, one may obtain the irrotationality of the peculiar-acceleration by time differentiating 
Eqs. (2.52-2.53). 

2.5 Remark 

Before we discuss the solutions of the aforementioned sets of equations we wish to say a few 
words on their usage in the current literature. 

Eqs. (2.29-2.31) are also reported in [38], however the irrotationality condition for the 
peculiar-velocity (and hence for the peculiar-acceleration as well) is flawed. Specifically, once 
the tensor T^ ■ • is nonzero, their expression for T^ n ' is false. T^ • (q) ^ means that 
the inhomogeneous deformation tensor Vl/j ,(#,77) consist also of an antisymmetric part, i.e., 
in the most general case T^ ' ^ and thus T^ + ' is not correct, with k £ No- 

The corrected irrotationality condition for the peculiar- velocity in [38] must read: 

s ijk [(1 + *n,«) S tj - *,j + Vfj] tf' M = , (i = 1, 2, 3) . (2.54) 

Eg. (2.54) is equivalent to Eq. (2.29), as can be seen by decomposing an arbitrary tensor 9ij 
into a symmetric and an antisymmetric part and evaluating both equations. Note that the 
subsequent perturbative calculation simplifies clearly with the usage of Eq. (2.29). 

Finally, we would like to stress again, that Lagrangian and Eulerian transverse mo- 
tions are not the same, since both frames are connected by a non-linear transformation (see 
Eq. 2.28). In reference [42] the Eulerian irrotationality condition was not solved, and as a 
wrong consequence they concluded that the Lagrangian transversality is zero at all orders, 
as long as the motion is irrotational in the Eulerian frame. In order to avoid confusion with 
this issue in the following, we shall refer to 'irrotational' only with respect to the Eulerian 
frame, whereas we reserve 'transverse' to the Lagrangian frame. 

3 From general initial conditions to 'Zel'dovich type' initial conditions 

The general solution in Lagrangian perturbation theory consists of solving the above set of 
equations with (3+1) initial condition functions. However, for an irrotational fluid we have 
to specify only two initial potentials, one for the peculiar- velocity u(q, to) = V q S± (q, to), and 
the other for the peculiar-acceleration g pec (q,to) = V q 52(q, to) (where the latter is linked to 
the density contrast via the Poisson equation). 

In [3] it has been shown that under the consideration of the most general ansatz for 
the integral curves of the full- velocity field /, one may express the solutions in terms of the 
general initial conditions setting. For example, in the full-system approach the first-order 
integral curves /' ' read in terms of the initial conditions: 

/« = aq + bi(t) u D {q, t ) t + b 2 {t) u R (q, t ) t + b 3 (t) g° cc {q, t ) t 2 , (3.1) 



where we have decomposed the peculiar-velocity in a curl-free part u and a divergence-free 
part u , i.e., u = u +u , and similar for £7 pec , but in this formulation the vortical part 
of g pec is zero, thus g pcc = g pcc - The explicit time coefficients bi(t) are not needed for this 
demonstration, but see Eq. (23a) in [3]. 

To avoid formally lengthy expressions for the solutions we impose a functional de- 
pendence on the two initial potentials S\ and S2 ■ The simplest classes of such functional 
constraints, that we call 'Zel'dovich-type' initial conditions, only prescribe one initial poten- 
tial [15], e.g. S\ = S2 = S. In the following, we describe two different settings, which are 
consistent within the Zel'dovich approximation and are thus appropriate for studying the 
inhomogeneities of the large-scale structure: the slaved initial conditions and the inertial 
initial conditions. 

Slaved initial conditions. We require the following parallelity condition between the 
initial peculiar-velocity u(q, to) and the initial peculiar-acceleration g (q,to) (note that, in 
general, parallelity allows for an arbitrary scalar proportionality function): 

u{q,t ) =g pcc {q,t )t . (3.2) 

From this it follows that u D (q,to) = g pcc (q,to)to and u R (q,to) = g pcc (q,to)to = 0. Thus, 
the integral curves (3.1) reduce to the expression 

f { z\=aq+[b 1 {t) + b z {t)} u D (q,t )t . (3.3) 

Since g pec (q,to) is non- vanishing, we have V q • g pcc {q,to) oc 5o initially, where the initial 
density contrast figures in the exact expression S = (1 + 8q)/J f — 1. 

In [15, 16, 43] it has been pointed out that the above restriction is appropriate for mod- 
elling the large-scale structure: gravitational instability supports the tendency for density 
perturbations to grow along potential flows — as long as one restricts the problem to irrota- 
tional flows. These so-called slaved initial conditions are thus a dynamical outcome of the 
Eulerian linear perturbation theory, where the growing mode supports the above parallelity 
condition. Consequently, this class of initial conditions is well-motivated, and under this 
initial assumption parallelity is then exactly preserved by the first-order Lagrangian solu- 
tion. Physically, this setting implies initial density- and velocity perturbations, which is in 
accordance with the prediction of standard inflationary theories. This initial condition is also 
commonly used in iV-body simulations. 

With the constraint Eq. (3.2) we only need to specify one potential initially, and set 
w(q,to) : = Vq5(q,to) and, thus, 

ffpec(<7,*o) = V q S(q,t )t . (3.4) 

According to the initial density-perturbation which is specified to be po = ~p Q (1 + e<5o) we have 
to relate the initial potential S(q, £0) to the potential <p ' (q) of the longitudinal perturbations: 
^ (q) — Vq4- '(<?)■ We shall do so by using Poisson's equation for Eq. (3.4). The slaved 
initial conditions shall be used in the following section 4. 
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Inertial initial conditions. Another class of initial conditions is to require that 

9 P ec(q,to) = 0, (3.5) 

which can be considered as a 'Zel'dovich-type' class as well. From this it follows that 
u R (q,to) = and gp ec (q, to) = 0. Then, the integral curves (3.1) become 

f { l\ =aq + b 1 (t)u D (q,t )t . (3.6) 

We call this class inertial initial conditions, since the initial acceleration is vanishing. In 
contrast to the slaved initial conditions this implies that there is no initially prescribed 
density perturbation, because Poisson's equation dictates a vanishing source term initially 
(see Eq. (4.3)). Instead we demand an initial homogeneous background density ~p(to), and 
the evolving density perturbations are produced solely by initial velocity perturbations [15]. 
As before we can then set initially the peculiar- velocity to u(q,to) := V q «S(g,io)- hi this 
case, mass conservation is given by the exact setting S = 1/J — l. 4 We shall use inertial 
initial conditions to solve the equations in the peculiar-system approach (see section 5). 

Although inertial initial conditions are of the 'Zel'dovich-type', it is difficult to justify 
their physical motivation. The main reason is because it is nowadays believed that the CMB 
fluctuations are of adiabatic nature, and the (non-relativistic) counter-part of adiabatic initial 
perturbations would rather correspond to slaved initial conditions. However, as we shall see 
in the upcoming sections, only decaying modes change, whereas the fastest growing modes 
coincide with the use of both initial conditions. 

In an expanding universe, the physical significance of the decaying modes is usually 
negligible, because they are dominated by the fastest growing mode. However, in the highly 
non-linear regime, decaying modes can affect the clustering via their coupling to growing 
modes [16]. (Also, since the equations and solutions are time-reversable, their application 
backward in time exchanges the role of the importance of growing and decaying modes.) 
Restricting the general initial conditions setting is mainly motivated by simplification. There 
is a price to pay for these restrictions when one is interested in the structural details at highly 
non-linear stages at small scales, but this regime lies in the realm of A-body simulations. 

4 Solutions of the perturbation equations for an EdS background: 
full-system approach 

In this section we solve the perturbation equations with the use of the slaved initial conditions, 
see Eq. (3.2). 

4.0 The zero-order solution 

To find the solution for the homogeneous and isotropic background, Eq. (2.41), one has to 
solve the Friedmann equation for the cosmic scale factor a(t). In the case of an EdS universe 
we have: 

/ t \ 2/3 
„«>=(-) . (4.1) 



4 With slaved initial conditions we may, however, assume an initial quasi-homogeneity, and hence approxi- 
mate mass conservation by 8 ~ 1/J F — 1, because of numerical smallness of the initial conditions. Alternatively, 
we may choose a different set of Lagrangian coordinates to render this form of the integral exact, but this 
will introduce formal complications. Zel'dovich's original model implicitly supposed slaved initial conditions, 
while the initial density contrast was approximated by zero in the density integral. 
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Using the Friedmann equation again we can write the prefactor AirGp^ in Poisson's equation 
as 2/(3tg). 

4.1 The first-order solution 

Using the above result and the longitudinal ansatz ^ l '{q) = Vq^W(qf), Eq. (2.42) reads: 

^ +2 H a ^ (i)= -4^'°- 

Comparing with Poisson's equation for the initial potential S(q, to), namely 

2 



(4.2) 



A q S(q,t )/t 



o 



3t 2 6 o = Vq-flf pcc (g,t ) 



(4.3) 



we can read off the connection between the spatial dependence and the initial conditions. 
The solution is then 



q\ + 2-qi 



^ ; A^iq)^^ =A q S(q,t )t 



aHl 



(4.4) 



The solution of the time-evolution under the restriction Eq. (3.2) can be found by the formal 
requirement qx(to) =0, and the condition that the coefficient functions of the peculiar- velocity 
and peculiar-acceleration are equal to 1 at t = to [16]: q\(to) = 1/to- The temporal solution 
is then 



qi = 2 ^ ~ a ] ' 

4.2 The second-order solution 

After separating the spatial part from the temporal part of Eq. (2.43) we have: 

[a 2 q 2 + 2aaq 2 ] = - [2qiq\a + aqj] ; fif = /4 . 
This linear ordinary differential equation can be solved and is found to be 



(4.5) 



(4.6) 



q-2 



3 q 6 n 8 1 In 

-a i + -a 2 -a+—a 1/2 
7 5 35 



(4.7) 



with the restrictions q2(to) = and 92(^0) = 0- The same restrictions are used for the higher- 
order solutions as well. 

4.3 The third-order solution 

Eq. (2.44) consists of two separate spatial parts with its temporal parts, the first of which 
describes a coupling between first- and second-order perturbations, the second is dependent 
on cubic first-order perturbations. As suggested in [15] the third-order displacement is split 
into two longitudinal parts, arising from these two specific dependencies and then solved 
piecewise. Thus, the solution of Eq. (2.44) consists in solving the two sets of equations, 
which we label with 3a and 36. Additionally, we have the occurrence of a transverse solution 
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at the third order from Eq. (2.46), which we label with 3c and separate it into spatial and 
temporal parts as well. The three sets of equations then read: 



936 + 2-g 3 6 
a 



q 3a + 2-q 3a 
a 



,.. 9? 



a 92 . 9i 

-^qm + qi — I — 92 



(I 



a 



■T 
Q3c 



a T 
7fe 



9i ■ 9i 

—q2 92 

a a 



(3a) 

Mi 



M? 6) 



& 



,r ] 



(3c) 



J £ ijkT^j 



%**, ( >S ■ 



(4.8) 



As mentioned earlier the purely longitudinal parts are fi^ = A q (jy- n > . The purely transverse 
part can be written in terms of a vector potential, i.e., T^ = Ekimd qi A m , where one degree 
of freedom may be cancelled by employing the Coulomb gauge. Using the same restrictions 
as before, we obtain for the temporal parts of the three above equations: 

„ _ (3\3 r 1 4 . 9_3 9 _2 i „ 16 .,-1/21 
93a - {2) L~3° + 7 a - 5° + a - T05 a J ' 



_ /'3\ 3 T10„4 66 „3, 14 „2 o„ , 16_l/2 , 16 _-l/21 
936 - UJ l2T a ~ 35° + T a ~ 2a + 35 a + T05° J ' 



(4.9) 



T 



(|) 3 [-ia 4 + fa3 + ia 2 -a + faV 2 -^ a -i/ 2 ]. 



4.4 The fourth-order solution 



The fourth-order part of the source equation, Eq. (2.45), consists of five separate spatial 
parts with their respective temporal parts. Analogous to the last section we label them with 
4a— 4e. The fourth-order part of the irrotationality condition, Eq. (2.47), consists of two types 
of spatial parts, the first type arises from longitudinal perturbations only, i.e., £ijk^\ j^il> 
whereas the second type results from a mixing term of longitudinal-transverse perturbations, 



ie £ --,* (1) T (3) 



Due to the fact that the third-order longitudinal perturbations include the 
solutions 3a and 36, we obtain two terms contributing to the first type, which we label with 
4/ and 4g. The second type occurs only once which we label with Ah. The eight sets of 
equations to be solved then read: 



94a + 2-94 a 

a 



946 + 2-946 

a 



,-9l + 2^92 



92 2 , 9l 

2a^ + a^ 1(?2 



a 




94c + 2 -94 C 


= -2 


a 


- 


q^d + 2-94d 


= -2 


a 





94e + 2-94 e 

a 



a ,9i . 91 .. 

"^9l93a H 93a H 93a 

a . 9i . 9i .. 

-^9l936 H 936 H 936 

a z a a 

T 1 91 T , 9l --T 

-29i9 3 + —9 3 + —9 3 



(4a) (2) 

; ^1 = M 2 



1 - eu m e ipq V , y\>v j * 



; Mi 



(4c) (1,3a) 



. (4d) _ (1,36) 

, /x 1 — fi 2 



. „(4e) _ ^(l) T (3c) 
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and 



■T 
34/ " 


d T 

" -94/ 


= 


9i . 

— 93a - 

a 


9i 

— 93a 

a 


e--iT {Af) 
> t vk- L k,j 


= ««*•&»•£?. 




a T 
a y 


= 


9i • 

— 936" 
a 


9i 

— 93b 

a 


. - T {4g) 


-««.»!>S>. 


■T 

94h" 


a T 


= 


"91 -r 
— 93c ~ 

a 6C 


9i t" 

— 93c 

a 6C 


' t -i-jk J -k,j 


-e ijk m l:j i lk . 



(4.10) 



Using the same restrictions as before, we obtain for the temporal parts of the eight sets of 
equations: 



94a 



946 



9ic 



94d 



9le 



T 

qi f 



T 
lA 9 



T 



(I) 4 [■ 



. «_ a 5 , 4 4 _ 258 3 , 12 2 _ 48 3/2 
539 " T 7 U 175 " ^ 5 " 49 " 



a +175V a +385 a 1225° J ' 



1/21 



(I) 4 [it* 5 " it* 4 + fa3 - fa* + ia^ + i a - l^ - ^ a -x/ 2] , 

(I) 4 [M« 5 " If « 4 + W« 3 " 4« 2 + 2a - #^ - §,a^ 

(I) 4 

(I) 4 

(I) 4 

f^ 4 rj_„5 116 _ 4 , 66 n 3 4 2 48 3/2 , , _32_ /- 16 --1/21 
UJ L2T a 105° +35° 5 a 35° + a + 105 V a 105° J 

(I) 4 hlV + f ° 4 " f° 3 + h 2 ~ f aV2 + f«- T^+ i a " V2 ] • 



["I* 5 + i« 4 " W« 3 + f « 2 " 4« + » V~a + j& 



(4.11) 



[- 



!„5 J_ IS _4 _ 9„3 _i_ 4^,2 _,_ 1 
6 U "^ 21 a 



-L 4 ^ 2 -^ 1 ^ 16 /77 _|_ 16 n-1/21 
+ 5 a +2 a_ 2TVa+TM a 



5 Solutions of the perturbation equations for an EdS background: 
peculiar-system approach 

As promised above, we solve the perturbation equations in this section with inertial initial 
conditions, see Eq. (3.5). 

5.1 The first-order solution 

The first-order part of the peculiar-source equation, Eq. (2.48), is 

D"{r,)^ q ^=a(r 1 )D(r 1 )^\ (5.1) 

with a = 6/tj 2 . We fix the constants of the temporal solution D{rf) such that the normal- 
isation coincides with the one obtained by linear SPT. Formally this can be achieved by 
the requirement that D(rjo) = 0, and that the coefficient function of the peculiar- velocity 
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is initially equal to —5: u(q,rjo) = l/a(r]o)D'(r]o)'V q 4>( 1 > 
then 



■■^■VJW. The full solution is 

Vo q ^ 



1 



-A q S(q,f] )r] Q 



D = a-a~ 3 / 2 ; A^«(g) = M « 
with the initial condition u(q,rjo) = 'V q S(q, t]q). 

5.2 The second- order solution 

Using Eq. (5.1) we can rewrite the second-order solution of Eq. (2.49) as 

/4 2) = /4 1} ; E " -aE= -aD 2 . 



Then we obtain 



E 



3 2 5 

~^ a + 7 a 

7 4 



20-V2 + H„-3/2 



1 



(5.2) 



(5.3) 



7 - 4- ; _ (5 ' 4) 

with the restrictions E(tjo) = and E'(j]q) = 0. The same restrictions are used for the 
higher-order solutions as well. 



5.3 The third-order solution 

Analogous to section Eq. (4.3) the third-order solution consists of two longitudinal (a and b) 
and one transverse solution (c). Using the lower-order results we can write the solutions of 
Eq. (2.50) and Eq. (2.52) as follows: 



F" - aF a = -2a D 3 



F' h ' - aF b = -2aD [E - D 1 



F" 



-aD' 






u (1) 

(1,2) 
M 2 



(5.5) 



.(3c) 



s^h® 



£ iJk T k,j fc "«* *l'J *J,fc 

For the three time-evolution functions we obtain under the same restrictions as above: 
F a = - la 3 + ga - 9aV2 + 25 a -3/2 _ 6a -2 + ^ a -9/2 ? 



A a -9/2 
22" 



+§a" 2 + fa" 3 
Fj = -la 3 + 9aV2 + «g> a -i/2 _ - a 



(5.6) 



3„-2 , J_ n -9/2 _ 125 
"i" 12" 6 



5.4 The fourth-order solution 



The solutions of the source equation, Eq. (2.51), and of the irrotationality condition, Eq. (2.53) 
can be written as follows: 



G" a -aG a = a [2D 2 E - E 2 



Gl - aGb 



a 



D 2 E 

2 



(4a) 

Mi 

,(46) 



(2) 
M 2 



(!) ,!,(!) ,T,(2) 



Ml — £ tkl£jmnV km V ln <Xj 



G" c - aG c = -2a [DF a - 2D 
G" d - aG d = -2a [DF b - 2D 2 E + 2D 



G e — aG e 



G 



rpll 



U g ■- 



-aD q 
2aD 4 

21 



-2aD 2 [E- D 
-aD [D 3 + F ( 



M 4c) 

Mi 4d) 

M 46) 



..(1,3a) 

P2 ) 

M 2 X ' 36) , 

*(l) T (3c) 



(5.7) 



£ *3k 1 k,j — e Vk *l,j *f,Aj ) 






(49) 






(4h) 
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We then obtain for the temporal parts of above sets of equations: 

r - --5Ln 4 4- 25« 3 - -15-n 2 - An 3 / 2 4- 2125,. _ 300-1/2 , 25-1/2 _ 51 _-l 
a — 539 42 112 ~ r * 231 49 ~r y u ^u, 

_250 --3/2 i 25-2 _ 25-3 _ 2„-7/2 , 25-9/2 3__-6 

231 " ' 8 " 49" 3 U ' 77" 112" ' 

r\ J^-n 4 - A„3 _ 6-3/2 , 2375 - _ 45-1/2 , 27--1 _ 125 --3/2 

^b — 154" 24 u 7" " r 924 " 28" " r 56" 66 " 

~"~28 a ~~ T2 a ~~ 77 a ' 28 a ' 

Gc = l| a 4 _ ^) a 2 + 4 a 3/2 + m a + HX) a -l/2 _ 27a -l + 1500 a -3/2 
+ f a -3_52 a -7/2 + 16 a - 6) 

r — _20 n 4 , 25-3 , 150-2 , 80-3/2 _ 4000- , 75-1/2 _ 1150- -1/2 , 75--1 
°" ~~ 33 u "i" 21 " "•" 77 " "r 21 231 ^ 7 231 ' 7 

+ |0 a -3/2 _ 50 a -2 _ 50 a -3 + g a -7/2 + 50 a -9/2 _ 20 fl -6 > (5 . 



G e =-ia 4 + 4a 3 / 2 - ^ a + Qa^ 1 - f^a" 3 / 2 + §a" 



7/2 _ 1 _-6 



21 



G J = -la 4 + 4a 3 / 2 + fP«- 1/2 - 36a" 1 + f 



n-T/2 _ J_ -6 _ 625 
" 11 " 14 ' 



nT — 5 _4 5 -3 20 „3/2 , 45„l/2 8125- -1/2 , 225 _-l , 45 _-2 
^9 _ 2T a ~~ I4 a ~~ T ° + T a ~~ ~23T a + TT a + 28 a 

_5 a -7/2 _ A a -9/2 , A a -6 , 625 
7" 21 " T 44" T 42 ' 

n T - 1 n 4 !8„3/2 , 125- 1125 -1/2 , 27--1 125--3/2 
Lr h -~i4 a _ T° + ~6 _a__ 28 _a + T a ~ ~6~ a 
, 375 --2 , 29 --7/2 _ J_„-6 , 625 

' 56 " " r 84" 24" "i 28 ' 

respectively. 

The above solutions will be used in the following section. For our purpose, i.e., setting 
up the relation to SPT, only the fastest growing modes will be relevant. 

6 Practical realisation and exact relationship between LPT and SPT 

In this section we explicitly link the perturbative series between LPT and SPT. Note that 
we still restrict to an EdS universe and to the Zel'dovich-type subclasses of the general 
Lagrangian solutions. As mentioned above, these subclasses consist of the slaved and the 
inertial initial conditions. 

Commonly, the dynamical quantities of the irrotational solutions in SPT are the density 
contrast and the peculiar-velocity divergence, so it is convenient to use in the following the 
peculiar-system approach only. In the latter we call \1/ the (gravitational induced) displace- 
ment field. 

6.1 The LPT series in Fourier space 

We seek a series solution of the form 



*(k,t) = Y,[* {n \k,t) + f (n) (k,t)} , (6.1 



n=l 
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where a tilde denotes Fourier quantities, and we explicitly neglect the perturbation parameter 
e. Eq. (6.1) is the complementary series of Eq. (2.35) in Fourier space for the displacement 

field. Then, using the solutions of section 5, it is possible to write the Fourier transform 

- (n) - (n) 

of the spatial nth-order perturbations \l/ (k) and T (k) times the fastest temporal parts 

oc D n (t) as 5 

#<■>(*,«) = -iZ J "( t )/ d3 ^ )3 d > (2 J )^g»( Pl ...„- fc ) 

xS<"'(p 1 ,...,P„)*< 1 '(Pl)---* (1) (P„), (6-2) 



and 



T W (M) = -U3"(t) / ^[^ (^M\ P ,..n - k) 



xS? ) (p lr .,pJP(p 1 )-PW 1 (6.3) 

where we have employed the shorthand notation P\... n = Pi + p 2 + • • • + P n , <!> and Sy 
are symmetrised vectors arising from longitudinal and transverse perturbations, respectively 
(see appendix A. 3). S^'(p) is the Fourier transform of the linear density contrast evaluated 
at to- We derive the above expressions in detail in appendix A. 

It is important to be conscious of the origin of the linear density contrasts on the RHS in 
Eqs. (6.2-6.3), since they include a specific limit. To see this we revisit the mass conservation: 

S(x,t) = jF -1, with x(q,t) = q + *(q,t). (6.4) 

Linearisation of this (for any displacement exact) Lagrangian expression implies for the spa- 
tial part 

6U(x)K-A q cf>U(q), (6.5) 

with \l/^ '(q) = Vq^' 1 ^^). In the ZA, the displacement field is very small, so it is appropriate 
to use the limit x ~ q, however neglecting the inherent non-linearity of the Lagrangian 
approach in this step. Thus, the quantities x and q in Eq. (6.5) are interchangeable, and this 
reads in Fourier space 

|W(p)«/<p)(p), (6.6) 

~ (n) ~ (n) 

at the linear order. The n linear density contrasts in \1/ and in T result from n times the 

~ (n) - (n) 

usage of Eq. (6.6). Stated in another way, * and T are affected by this approximation 
at any order, but strictly speaking, only because the above linear approximation is used n 
times. We call this setting the initial position limit (IPL). 6 

In [5, 44] the same setting was used implicitly in order to derive the longitudinal solution 
Eq. (6.2) up to the third order in LPT. 



5 D(t) is the fastest growing mode solution of the linear growth function in the peculiar-system approach, 
see Eq. (5.2). For an EdS universe it is D(t) = a(t). 

6 Given an arbitrary tensor function T(x,t) with the trajectory F = x(t) = q + *P(q,i). Then the IPL is 
Tipi J (x,t) := lim*_).o T(x, t) = T(q,t), which is a first-order approximation to T(F^ 1 (x,t),t). 
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6.2 Exact relationship between the density contrast and the displacement field 

In SPT, the series of the Fourier transform of the density contrast can be formally written 



as 



5(k,t) = J2~^ n) (k,t). (6.7) 



n=l 



Solving the Eulerian equations of motion order by order in SPT, the resulting nth-order 
density contrast in Fourier space reads: 

#»>(M) = D n (t) f d3pi {2n)3 d y n (2vr) 3 5g ) (p 1 ... n - k) 

xF( s )(p 1) ...,p n )P(p 1 )-^W, (6.8) 

(s) 

where the Fn 's are the symmetrised SPT kernels given by the well-known recursion relation 
[8, 45], and are explicitly given in [33] for n < 4. 

In order to find the link between the above series and the displacement field, we use 
again mass conservation, i.e., Eq. (6.4). In Fourier space this is 

S(k, t)= J d 3 x e ik - x 5(x, t)= I d 3 q e lkq {e lk *^ t] - l) . (6.9) 

To get the last expression we used d 3 x = \dx/dq\ d?q, with \dx/dq\ = J F . Taylor expanding 
the above equation and explicitly using the LPT series, i.e., Eqs. (6.1-6.3), we can sum up 
all specific contributions and associate them to their respective nth-order density contrasts. 
The result is then 

~ § (n) {kit)=D nJ d ^ )3 d > ( 27r) 3,g) {pin _ k) 

xXW(fe;p 1 ,...,p n )|( 1 )(p 1 )...|( 1 )(p n ), (6.10) 

(s) 

where we have defined the symmetrised kernels X\ s, which can be found in appendix B for 
n < 4. The Dirac-delta in Eq. (6.10) fixes k = P\... n and, as a result, 

Xi s \k; Pl ,...,p n ) =FW(p 1 ,...,pJ, (6.11) 

k =Pl ■ n 

valid at least up to the fourth-order in perturbation theory. As an important conclusion, 
Eq. (6.9) leads to identical expansions either in SPT or LPT. 

6.3 Exact relationship between the (peculiar-)velocity divergence and the dis- 
placement field 

Similar considerations can be made for the peculiar-velocity divergence. In Eulerian space 
the peculiar-velocity is u = dx/dr, where dr = dt/a is the conformal time, and as before 
x = q + \I/. For an EdS universe the time dependence of the fastest growing mode is oc a n to 
nth order in perturbation theory, therefore d/drE^)^^ '(q) = 2T-LE{t)^^ {q), and similar 
for higher orders, where T~L is the conformal Hubble parameter. Furthermore, the divergence 



7 The expression 8(k,t) = Jd 3 gexp{ifc-q + ife- *}[1 — J F ] is identical to Eq. (6.9), and leads to a somewhat 
different expansion but equivalent results. 
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of the peculiar-velocity is defined by V^ • u(x,t) = 0(x,r). Up to the fourth order in LPT 
the peculiar-velocity divergence is then 

6(x, t) =7iDV x - * (1) (g) + 2UE V x ■ * (2) (q) + ZUF V x ■ * (3) (q) 

+mGV x -V^(q). (6.12) 

To proceed it is useful to transform the divergences into Lagrangian coordinates, i.e., d/dx = 
\dx/dq\~ 1 d/dq. Then, Eq. (6.12) translates to 

^^1 = J F [J^ (D^(q) + 2B»g(g) + 3F^(q) + 4G9${qj) , (6.13) 

where [JE]~ l = l/(2J F ) £um £ jpqJpiJ£mi with J F = 5ij + ^ij- As before commas denote 
derivatives with respect to Lagrangian coordinates q. The Fourier transform of 8(x,r)/7i is 

U(t) ~ J q H(t) ' [b j 

where we have used d 3 x = \dx/dq\ d 3 q, and J f 9(x,t)/'H(t) is given by Eq. (6.13). We seek 
a perturbative solution of the form 

oo 

e(k,r) = J20 {n) (k,r). (6.15) 

n=l 

Taylor expanding Eq. (6.14) and using the LPT results, i.e., Eqs. (6.2, 6.3), we obtain for the 
nth-order peculiar-velocity divergence: 

*»(fc,r) = -Hirf d3 ^ )3 f P " (2nM\ Pl ... n - k) 

x yW(fe;Pi,--- 5 p»)* (1) (Pi)"-« (1) (Pn), (6-i6) 

where the symmetrised kernels Yn up to the fourth order are given in appendix B. Then 
similar to the above we find 

Y^(k; Pl ,..., Pn ) = G$( Pl ,..., Pn ), (6.17) 

fc=Pi...„ 

(s) 

with Gn given by SPT recursion relation. Thus, the prediction of the peculiar-velocity 
divergence agrees with its counterpart in SPT, which is: 

6^(k,t) = D»(t) J ^^ (2vr) 3 4 3 V... n " fc) 

xGW(p ll ...,pjl 1 )(p 1 )...|«(p„). (6.18) 

The SPT kernels Gn up to the fourth order are given in [33] as well. 
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7 Discussion and conclusions 

We have calculated solutions of the Lagrangian perturbation theory (LPT) for an irrotational 
fluid up to the fourth order. The derivation is shown in two separate approaches, and 
we call them the full-system- and the peculiar-system approach. They are solved for two 
sets of 'Zel'dovich-type' initial conditions: for the full-system approach we use the slaved 
initial conditions, where the initial parallelity between the peculiar-velocity and the peculiar- 
acceleration is assumed; for the peculiar-system approach we use the inertial initial conditions, 
where the initial peculiar-acceleration is set to zero. Both scenarios are appropriate for 
studying the large-scale structure. 

Solutions up to the fourth order have been found for the first time in [22] (in the full- 
system approach), whereas one of the current authors independently derived the solutions in 
the peculiar-system approach. 

In section 6 we transform the nth-order solution into Fourier space, see Eqs. (6.2-6.3). 
These solutions contain only the fastest growing mode and consist of integrals over n linear 
density contrasts. It is important to note that these results are derived in the so-called 
initial position limit (IPL), in which it is assumed that the aforementioned linear density 
contrast 5^\x) is evaluated in the vicinity of the initial Lagrangian position, i.e., in the 
limit x ~ q. Stated in another way, since the displacement field is small in the Zel'dovich 
approximation, we can evaluate the results in this limit, and the resulting Poisson equation 
for the displacement potential assumes therefore the linearised form Aq^ 1 ' {q) ~ — 5^'(q). 
Working in the IPL requires a small displacement, the starting assumption of the LPT 
series, but it nevertheless implies that we neglect the inherent non-linearity of the Lagrangian 
approach in this step. Strictly speaking, only the linear displacement potential and the linear 
displacement field are affected by this approximation, but due to the reason that the nth order 
displacement field is dependent on n displacement potentials, this approximation carries over 
to any order. This procedure is implicitly assumed in the current literature when working in 
Fourier space. 

Then, we express the Eulerian dynamical variables, i.e., the density contrast and the 
peculiar- velocity divergence in terms of the gravitational induced displacement field \l/: the 
predictions in Fourier space yield identical expressions for the LPT and SPT series at each 
order in perturbation theory in the IPL. However, this does not imply, that the convergence of 
the series in SPT and LPT is identical, rather the above subclasses of the general Lagrangian 
solution can be used to mimic the SPT series, as long as one restricts to the IPL. Additionally, 
the IPL suggests that there is no proper distinction between the Fourier transform on Eulerian 
space and on Lagrangian space. Relaxing the IPL would thus result in SPT and LPT 
solutions in separated Fourier spaces, whereas these spaces are connected via a non-linear 
transformation. This scenario may lead to an even better approximation of the Eulerian 
variables in terms of LPT. However, this issue is beyond the scope of this paper and we leave 
it as a future project. 
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A Preparing the solutions for practical realisation 

Our goal in this appendix is to derive Eqs. (6.2-6.3). To do so, we first have to transform 
our spatial nth-order solutions, V q • *f?( n '(q) = A q (j)^ n '(q) and V 9 x T^ n \ into Fourier space. 

~ (n) - (n) 

Then, in appendix A. 2 and A. 3, we obtain the full nth-order solutions \J/ (k, t) and T (k,t) 
by multiplying the fastest growing mode solution with its spatial nth-order solution. 

A.l Fourier analysis 

In this appendix we prepare the Lagrangian solutions for their use in Fourier space. Let us 
start with our Fourier convention, 

A(k)= fd 3 xe ikx A{x), A(x) = i -^ e~ [kx A(k) . (A.l) 

As mentioned earlier, the nth-order longitudinal perturbation ^^ n '(q) can be written in 
terms of a potential Vq</>' n )(q). In Fourier space this means 

-ipi# n >(p) = *! n) (p). (A.2) 

We shall first focus on the net longitudinal contributions. The spatial part of the nth-order 
perturbation, e.g. Eqs. (5.3, 5.5, 5.7), yields: 

XK»(Pl, ■ ■ ■ ,P„) 4> m (Pl) ■ ■ ■ i m (Pn) ■ (A.3) 

In above equation we employed the shorthand notation Pi... n = Pi + p 2 + ■ ■ ■ + p n , and the 
kernels kS 71 ' are given by: 

^(Pl,P 2 ) = g [^2 - (Pi • P2) 2 ] , 

K3a (Pi , P2 , P3 ) = ~ g £ ikl PliP2kP3l £jmn PljP2mP3n , 

I n / x K 2(Pl,P 23 ) 

K3b{Pl,P2,P3) = - K 2{P2,P3) 2 ' 

P23 

«4a(Pl,P 2 ,P3>P4) = K 2(Pl,P 3 )«2(P 2 ,P4)^-^V^^' 

P13P24 

K46(Pl»P2»P3>P4) = ~6«3a(Pl,P2»P34) 2 ^ ' ( A - 4 ) 

^34 

«4c(Pl,P 2 ,P3.P4) = -^3a(P 2 ,P 3 ,P 4 ) — o 234 . 

^234 

/ \ / N K 2(Pl,P 2 34) 

«4d(Pi,P 2 ,P3>P4) = - K 3&(P 2 >P3;P4) 2 ' 

£*234 

K4e(Pl,p2»P3.P4) = ~ [(Pi ' P2)(P234 • P34) ~ (Pi ' P34)(P234 ' P2)] 

X (Pl-P234) P 2 P 2 4 k 2(P 3 ,P 4 )- 
^234^34 
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It is worthwile to make the derivation of K^ e explicit. In order to write an expression for <f>^ e \ 
we first note that the Fourier transform of the transverse part (last expression in Eq. (5.5)) 
can be written as 

, / „ ^ _ / fpM* (2l) 3 4 3, (p _ (p^p*) x 

V Ji J (2vr) 9 p| 3 

X«2(p 2 ,P3)^ (1) (Pl)^ (1) (P 2 )^ (1) (P3)- (A.5) 

Since T' ' is purely transverse, we can write it as a curl of a vector potential, i.e., T- = 

/■q\ 

ZjlmAmj. Using this, Eq. (A.5) reads: 

/if > - „ ( P . a"') = - / ^y* (2.)=4 3 >(p 123 - P ) <^p> 

x(Pi XRj), K2(p 2 ,P 3 )0 <1) (Pi)* ll, (P2)* ll, (P 3 )- (A.6) 

The second term on the LHS may be removed by an appropriate choice of gauge, however, 
this term will vanish anyway. Using Eq. (5.7), we can write: 

^ (4£) (P) = ^ / ^^ &*)* 4 3) (Pl2 " P)(Pl • P 2 )Pl^7mW (1) (Pl) (A-7) 

f 1 f d 3 p 3 d 3 P4 d 3 p 5 3 (3) x (P 3 -P 45 ) 

x(P3XP45) m ^(P4,P5)^ (1) (P3)^ (1) (P4)^ (1) (P5)}+^(P 2 -^ (3) (P 2 )) 

J P 2 

which yields K 4e . 

Now we proceed with the derivation of the nth-order transverse perturbation. The 
transversality requirement of T^ n > in Fourier space reads: 

ff i (jp) = -ie jlmPl At ) (jp). (A.8) 

Using this and the transverse results in Eqs. (5.5,5.7), we can write 

T (P) = -^J (27r) 3n ( 27F ) ^(Pl-n-P) 

x^ n (p 1 ,...,pj0( 1 )(p 1 )...0( 1 )(p n ), (A.9) 

where we have defined: 

<*>3c(Pl,P2,P 3 ) = [Pl(Pl 2 3-P23)-P 2 3(Pl23-Pl)] (Pi " P 2 s) ~ ^'^ » 

P23 

<*4f(Pl,P2,P3iP4) = [Pi (Pl 2 34 • P 2 34) ~ P 2 34 (Pl 2 34 "Pi)] (Pi " P 2 34) "^ ~\^ 



P 2 34 

U4g(Pl,P2,P3,Pi) = [Pi (Pl 2 34 ' P 2 34) ~ P 2 34 (Pl 2 34 " Pi)] (Pi " P 2 34) h3bP2 f 3,P4 , 

P234 

TJ ■ (ij ( Ti TJ T) 1 

W4A(Pl,P2.P3.P4) = [Pl(Pl 2 34-P 2 34)-P 2 34(Pl 2 34-Pl)]^ 3c_2 1 _^^^_ _ ^ jq) 

^ 2 34 
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It is important to note that the u;'s satisfy the condition: 

Pl2-n- u >n(Pl,---,Pn) = °- (A- 11 ) 

In the following it is also useful to symmetrise the above kernels. The symmetrisation pro- 
cedure is [45]: 



4 S) (Pli ■ ■ ■ ,Pn) = ^_Y1 K "(Pi(i)> • • -'Pliny) , 

' ies„ 

^n\Pl,---,Pn) = —^2^>n(Pi(l),---,Pi(n))- (A. 12) 



n! 

ies n 



A. 2 Time-evolution of the fastest growing mode for an EdS universe 

For practical use (that will become clear below), we write the time-evolution of the fastest 
growing mode in terms of the linear growth function D{t): 



D = 


= a 


n — 51 r>4 
Lxa ~~ 539^ 


G l = 


-±L> 4 


E = 


7 ±j 


G b = mD 4 


G T g = 


±D 4 
21-^ 


F a = 


- -in 3 


r< 14 n4 


r T - 
u h - 


14 ^ 


F b = 


21"^ 


C 20 n 4 

Gd - -33^ 






F c = 


= -±L> 3 
7 u 


G e = -^L> 4 







(A.13) 



For general cosmologies, VL m / 1, SI a / 0, one can also find approximate solutions of the 
time-evolution, i.e., for Eqs. (5.4, 5.6, 5.8). For the first-order up to the third-order solutions 
this is often performed by employing fitting factors. However, the impact of Q m / 1 on the 
perturbative kernels is very little and can often be neglected, and we assume that this is also 
the case for the fourth-order solution. 



A. 3 Perturbative displacement fields 

By using Poisson's equation, the linear displacement potentials in Eqs. (A. 3, A. 9) can be 
linked to (the spatial part of) the linear density contrast: p 2 (j)^ 1 >{p) ~ ^^'(p). This is the 
initial position limit. Then, all we have to do is combine the results from the last two sections, 

~ ( 2 ) ~ (2) 

e.g., the Fourier transform of the second-order displacement fields is \]/ (k, t) = E(t)*& (k). 
Using Eq. (A. 2) we finally obtain the (symmetrised) result: 

*<">(M) = -UP® J ^^ (2nM\ Pl ... n - k) 

xSW( Pl ,...,p n )~S ( - l \p 1 )---~6U(p n ). (A.14) 
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S^ n ' are the longitudinal perturbative vectors, themselves dependent on the symmetrised 
nth-order kernels k„ (see Eqs. (A. 4)): 

pi 



S {2) ( Pl ,P2) 



o (a) 

6 Pl2 K 2 

T 2 2 2' 
' Pi 2 P1P2 



S( 3 )( Pl ,p 2 ,p 3 )= 2 Pl f 2 2 
P123P1P2P3 



P1234 



3 



- K 3a 21 ^ 



(A.15) 



S (4) (Pi, 



P2,P3,P4) 



Pl234PlP2P3Pl 



51 ( s ) _ J^_ (s) 

539 * 4a 154 Kib 

_ 14 (s) 20 (s) 1 w 
33 4c + 33 4rf + ll K4e 



where we have suppressed the dependences of K n = K n (p 1; . . . ,p n )- 

Computing the transverse part T (k,t) is then completely straightforward. From 
Eq. (A. 9) and Eq. (A. 13) we obtain 



f {n) (k,t) = -iD n {t) 



d 3 Pl • • • d 3 p n < 
(2vr) 3n 



(2vr) 3 ^( Pl ... n -fc) 



4 n) ( Pl , . . . ,p n ) ~5^( Pl ) ■ ■ • *«(pj , (A.16) 



n), 



and the only non- vanishing S^, 's are 



? (3) 



1 u 



3c 



^T (Pl)P2)P3) ~~ v~2 2 2 2' 
7 PI 23 PW3 

'S't (Pl'P2>P3'P4) = ~ 2 2 2 2 

Pl234HP»4 



6 4 / " 21^ + 14^ 



(A.17) 



up to the fourth order. The oj's are given in Eqs. (A. 10). It is important to note that 
P\2-m ' &t (Pi) • • • >Pn) omv vanishes if m = n, due to transverseness. 



B Kernels for the relationship between SPT and LPT 

The symmetrised kernels for Eq. (6.10) are 

x[ s \k; Pl ) = k-SW( Pl ), 

X { 2 s \k; Pl ,p 2 ) = k ■ S (2) (Pi,P 2 ) + \k ■ S (1) (Pi) k ■ sW(p 2 ) , 
X^ s \k; Pl ,p 2 ,p 3 ) = k • S {3) (p 1 ,p 2 ,p 3 ) + S$\p 1 ,p 2 ,p 3 
+±k-SW( Pl )k-SU(p 2 )k-sU(p 3 ) 
+ - {fc • S (1) (Px)fc • 5( 2) (p 2 'P 3 ) + two perms.} , 



(B.l) 
(B.2) 



(B.3) 
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X { 4 s '(k;p 1 ,p 2 ,p 3 ,p 4 ) = k- S (4) (p 1 ,p 2 ,p 3 ,p i ) + S { T '{p 1 ,p 2 ,p 3 ,p 4 ) 
+±k ■ SW( P1 ) k ■ S^(p 2 ) k ■ 5«(p 3 ) k ■ S«(p 4 ) 

+ 3 | 2 k ' S(2) ^ p ^ k ' s(2) (P3,Pi) +two perms. I 



1 



+ j & • 5W(po fe • [s (3) (p 2 ,P3,p 4 ) + 4 ) (p 2 ,p 3 ,p 4 : 



+ three perms. 



+ ^ j 2^ ' SW{Pl) k ' S(1)(P2) k ' 5(2)(p 3'P4) + five perms. 



and for Eq. (6.16) 




Y} s \k; Pl ) = k.S^( Pl ), 


Y 2 {s \k; Pl ,p 2 ) = ~ 


\ (Pi'Pi) 2 
11 



+ lk.S^( Pl ) + lk.S^(p 2 ). 



Y 3 {s \k; Pl ,p 2 ,p 3 ) = Uk- S^( Pl ,p 2 ) + h ■ S«( Pl ) k ■ S«(p 2 ) 
-0P12 • S^ 2) (p 1 ,p 2 )k- S W {p 3 )+ two perms. I 



+ 3p 123 -S( 3 )- 3-5-5-5 + 



w (3o, 8 ) lg K (3b,a) 



(B.4) 

(B.5) 
(B.6) 



Y} S \k;p 1 ,p 2 ,p 3 ,p 4 ) 



P1P2P3 7 P1P2P3 

36 n^ a ^ 6 k( 46 ' s ) 8 k( 4c ' s ) 80 k^^ 



(B.7) 



+ 



49p{p$p$p( ipipippi ^p{p 2 p 3 pi ^Pipipwi 



'-4^ + *™ ■ ^ + i| - » ■ S «( PI )^^^ 



7P1P2P3P4 



p\p\p\ 



+h • ^'(p /^'^ + -U • sW( Pl ) fc • S«(p 2 ) fc • sW(ft) 
• P2P3P4 ° 



+fc 



+ 6 



s (3) (pi,p 2 ,p 3 ) + 4 3) (pi>p2,p 3 ) 



+ three perms. 



6^7 [ P3P4 J 7 L P3P4 

+ five perms. > + — < k ■ S^ '(Pi) k ■ S^ '(p 2 ,p 3 ) + eleven perms. > . 



(B.8) 



The validity of Eq. (6.10) and Eq. (6.16) with the kernels given in appendix B has been 
checked with Mathematica, and the respective code can be obtained upon request: 
rampfQphysik. rwth-aachen. de. 
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n 

t 

Qi 

V 

«(l) 

P(r,t) 

v(r,t) 

9(r,t) 

X i 

u{x,rj) 
9 pec (x,r)) 

m 

8(x,rj) 

f(q,t) 

F(q,n) 
J 

J F 
Po 

p(q,t) 
Pij 

*(g,r?) 

*(")(g) 

c 
9n(*) 

D(ij),E(rj), 

( n )/ \ 

WW) 

5(q,to) 

5 (n) 

o(") 



Eulerian coordinate (i = 1, 2, 3) Eq. (2.2) 

cosmic time (d/di = ') Eq. (2.2) 

Lagrangian coordinate (i = 1,2,3) Eq. (2.12) 

superconformal time (d/drj = ' ) Eq. (2.6) 

prefactor in the peculiar Poisson equation Eq. (2.6) 

density field Eq. (2.2) 

velocity field Eq. (2.3) 

acceleration field Eq. (2.3) 

comoving distance (i = 1,2,3) Eq. (2.6) 

peculiar velocity Eq. (2.10) 

rescaled peculiar acceleration field Eq. (2.24) 

mean density Eq. (2.5) 

density contrast Eq. (2.1) 

trajectory of the fluid particle Eq. (2.20) 

peculiar trajectory of the fluid particle Eq. (2.27) 

Jacobian in the full-system approach Eq. (2.13) 

Jacobian in the peculiar-system approach Eq. (2.26) 

initial density at q and to Eq. (2.16) 

perturbation in the full-system approach Eq. (2.32) 

cofactor element of the perturbation p Eq. (2.23) 

perturbation in the peculiar-system approach Eq. (2.35) 
(also called displacement field) 

cofactor element of the perturbation SI/ Eq. (2.31) 

longitudinal perturbation to the rath order (= Vg<p n )) Eq. (2.32) 

transverse perturbation to the nth order (= V q x A^ n ') Eq. (2.32) 

small and dimensionless perturbation parameter Eq. (2.32) 

growth function of the nth order in the perturbation Eq. (2.32) 
(full-system approach) 

linear growth function, second order growth function, • • • Eq. (2.35) 

(peculiar-system approach) 

ath scalar to the nth order Eq. (2.37) 

initial potential Eq. (3.4) 
perturbative vector arising from longitudinal perturbations Eq. (A. 15) 
perturbative vectors arising from transverse perturbations Eq. (A. 15) 

Table 1. Used notation. 



C Used notation 

In Tab. 1 we give an overview of the notation used in this work. 
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